/* libwatopt++ - Optimization Toolkit
 *
 * Copyright (c) 2007 Christopher Alexander Watford
 * <christopher.watford@gmail.com>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to
 * deal in the Software without restriction, including without limitation the
 * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
 * sell copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
 * IN THE SOFTWARE.
 *
 * $Id: kohonen_network.cpp 2 2007-05-06 18:01:15Z christopher.watford $
 */

#include <cstdlib>
#include <cmath>
#include <ctime>

#include <kohonen_network.hpp>

#include "utils.hpp"

namespace watopt
{
	double
	KohonenNetwork::Train(const double *inputs)
	{
		// find winner
		neuronid_t min_id = m_Neighborhood->First();
		double min_S = (*m_Neighborhood)(min_id, inputs);
		for(neuronid_t nn = m_Neighborhood->Next(min_id);
			nn < m_Neighborhood->Last();
			nn = m_Neighborhood->Next(nn)
			)
		{
			double ss = (*m_Neighborhood)(nn, inputs);
			if(ss < min_S)
			{
				min_S = ss;
				min_id = nn;
			}
		}

		// update winner and neighbors
      double tsse = 0.0;
		for(neuronid_t nn = m_Neighborhood->First();
			nn < m_Neighborhood->Last();
			nn = m_Neighborhood->Next(nn)
			)
		{
			tsse += m_Neighborhood->Update(nn, min_id, inputs, m_Rate);
		}

      return tsse;
	}

	neuronid_t
	KohonenNetwork::Categorize(const double *inputs)
	{
		// find winner
		neuronid_t min_id = m_Neighborhood->First();
		double min_S = (*m_Neighborhood)(min_id, inputs);
		for(neuronid_t nn = m_Neighborhood->Next(min_id);
			nn < m_Neighborhood->Last();
			nn = m_Neighborhood->Next(nn)
			)
		{
			double ss = (*m_Neighborhood)(nn, inputs);
			if(ss < min_S)
			{
				min_S = ss;
				min_id = nn;
			}
		}

      return min_id;
	}

	double
	OneDNeighborhood::Update(const neuronid_t nn, const neuronid_t winner,
		const double* inputs, double eta)
	{
		if(NULL == inputs)
			return 0.0;
		if(nn >= m_Neurons)
			return 0.0;
		if(winner >= m_Neurons)
			return 0.0;

		// not an error, but don't update non-neighbors
		if(!IsNeighbor(winner, nn))
			return 0.0;

      double tsse = 0.0;
		double* neuron = &m_Neighborhood[nn*m_Inputs];
		double weighting = Weighting(winner, nn)*eta;
      for(size_t ii = 0; ii < m_Inputs; ++ii)
      {
         double delta = weighting * (inputs[ii]-neuron[ii]);
			neuron[ii] += delta;
         tsse += delta*delta;
      }

      return tsse;
	}

	bool
	OneDNeighborhood::IsNeighbor(const neuronid_t n1, const neuronid_t n2) const
	{
		return (n1 > n2 ? (n1-n2 < m_HoodSize) : (n2-n1 < m_HoodSize));
	}

	double
	OneDNeighborhood::Weighting(const neuronid_t winner, const neuronid_t nn) const
	{
		if(winner == nn)
			return 1.0;

		double dist = (double)(winner > nn ? winner-nn : nn-winner);
		return ((double)m_HoodSize - dist)/(double)m_HoodSize;
	}

	double
	OneDNeighborhood::operator ()(const neuronid_t nn, const double *inputs) const
	{
		if(NULL == inputs)
			return -1.0;
		if(nn >= m_Neurons)
			return -1.0;

		double* neuron = &m_Neighborhood[nn*m_Inputs];
		double ss = 0.0;
		for(size_t ii = 0; ii < m_Inputs; ++ii)
			ss += (neuron[ii]-inputs[ii])*(neuron[ii]-inputs[ii]);

		return std::sqrt(ss);
	}

	void
	OneDNeighborhood::Get(const neuronid_t nn, double* weights) const
	{
		if(NULL == weights)
			return;
		if(nn >= m_Neurons)
			return;
		
		memcpy(weights, m_Neighborhood + (nn*m_Inputs), m_Inputs*sizeof(double)); 
	}

	void
	OneDNeighborhood::Initialize(void)
	{
		std::srand((unsigned int)std::time(NULL));
		
		for(size_t nn = 0; nn < m_Neurons; ++nn)
			for(size_t ii = 0; ii < m_Inputs; ++ii)
				this->m_Neighborhood[nn*m_Inputs + ii] = rand01();
	}
}

